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Abstract 

We present the derivation of a new unidirectional model for unsteady 
mixed flows in non uniform closed water pipes. We introduce a local 
reference frame to take into account the local perturbation caused by the 
changes of section and slope. Then an asymptotic analysis is performed to 
obtain a model for free surface flows and another one for pressurized flows. 
By coupling these models through the transition points by the use of a 
common set of variables and a suitable pressure law, we obtain a simple 
formulation called PFS-model close to the shallow water equations with 
source terms. It takes into account the changes of section and the slope 
variation in a continuous way through transition points. Transition point 
between the two types of flows is treated as a free boundary associated 
to a discontinuity of the gradient of pressure. The numerical simulation 
is performed by making use of a Roe-like finite volume scheme that we 
adapted to take into account geometrical source terms in the convection 
matrix. Finally some numerical tests are presented. 

Keywords : Shallow water, mixed flows, free surface flows, pressurized flows, 
curvilinear transformation, asymptotic analysis, VFRoe scheme, well-balanced 
finite volume scheme, hyperbolic system with source terms. 



1 Introduction 

The presented work takes place in a more general framework: the modelling 
of unsteady mixed flows in any kind of closed pipe taking into account the 
cavitation problem and air entrapment. We are interested in flows occurring 
in closed pipes with non uniform sections, where some parts of the flow can 
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be free surface (it means that only a part of the pipe is filled) and other parts 
are pressurized (it means that the pipe is full). The transition phenomenon 
between the two types of flows occurs in many situations such as storm sewers, 
waste or supply pipes in hydroelectric installations. It can be induced by sudden 
changes in the boundary conditions as failure pumping. During this process, 
the pressure can reach severe values and may cause damages. The simulation 
of such a phenomenon is thus a major challenge and a great amount of works 
was devoted to it these last years (see [I3])[Il])IlS],[in]) for instance). 

The classical shallow water equations are commonly used to describe free 
surface flows in open channels. They are also used in the study of mixed flows 
using the Preissman slot artefact (see for example [T31[5S]). However, this tech- 
nic does not take into account the depressurisation phenomenon which occurs 
during a water hammer except in recent works [2TJ [20l [22] . On the other hand 
the AUievi equations, commonly used to describe pressurized flows, are written 
in a non-conservative form which is not well adapted to a natural coupling with 
the shallow water equations. 

A model for the unsteady mixed water flows in closed pipes and a flnite 
volume discretisation have been previously studied by two of the authors [7] 
and a kinetic formulation has been proposed in [S]. We propose here the PFS- 
model which tends to extend naturally the work in [7] in the case of a closed pipe 
with non uniform section. For the sake of simplicity, we do not deal with the 
deformation of the domain induced by the change of pressure. We will consider 
only an infinitely rigid pipe. 

The paper is organized as follows. Section 2 is devoted to the derivation 
of the free surface model from the 3D incompressible Euler equations which 
are written in a suitable local reference frame (following O |4]) in order to 
take into account the local effects produced by the changes of section and the 
slope variation. The construction of the free surface model is done by a formal 
asymptotic analysis. Seeking for an approximation at first order gives the model 
called FS-model. In Section [31 we adapt the derivation of the FS-model to 
derive the pressurized model, called P-model, from the 3D compressible Euler 
equations. Writing the source terms of these two models, P and FS-model, 
into a unified form and using the same couple of conservative unknowns as in 
[8], we propose in Section [4| a model for mixed fiows, that we call PFS-model 
. We state some mathematical properties of this model. Section [5] is devoted 
to the extension of the VFRoe scheme described in O [161 [I] that was used for 
the case of uniform pipes. In Section [6l we show how to construct a convection 
matrix in order to get an exactly well-balanced scheme. Several numerical tests 
are presented in Section [71 
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Notations concerning geometrical variables 

• (0, i, j, k): cartesian reference frame 

• aj(x, 0, 6(x)): parametrization in the reference frame (0,i,j,k:) of the plane curve C which 
corresponds to the main flow axis 

• (T,N,B): Serret-Frenet reference frame attached to C with T the tangent vector, N the 
normal vector and B the binormal vector 

• X, y, Z: local variable in the Serret Frenet reference frame with X the curvilinear abscissa, 
Y the width of pipe, Z the B-coordinatc of any particle. 

• (t{X, Z) ^ /3(X, Z) - q(X, Z): width of the pipe at altitude Z with /3CX, Z) (resp. q(X, Z)) 
is the Y-coordinate of right (resp. left) boundary point at altitude Z 

• e{X): angle (i, T) 

• S{X): cross-section area 

• R{X): radius of the cross-scciion S{X) 

• J^wb^ outward normal vector to the wet part of the pipe 

• n: outward normal vector at the boundary point m in the Sl-plane defined below 

Notations concerning the free surface (FS) part 

• A: wet area 

• Q: discharge 

• Q(f, X): free surface cross section 

• H{t, X): physical water height 

• h{t, X): Z-coordinate of the water level, a-{X, h{t, X)) = T{A) : width of the free surface 

• rifs: outward B-normal vector to the free surface 

• pn: density of the water at atmospheric pressure po 

Notations concerning the pressurized part 

• 0,{X): pressurized cross section 

• p{t, X): density of the water 

• /3: water compressibility coefficient 

• c = —7^=: sonic speed 

\/0 PO 

• A = -^S: FS equivalent wet area 

• Q: FS equivalent discharge 

Notations concerning the PFS model 

• S: the physical wet area: S = A if the state is free surface, S otherwise 

• Ti: the Z coordinate of the water level: 7i = h ii the state is free surface, R otherwise 

Other notations 

• Bold characters are used for vectors except for S 
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2 Formal derivation of the FS-model for free 
surface flows 



The classical shallow water equations are used to describe physical situations like 
rivers, coastal domains, oceans and sedimentation problems. These equations 
are obtained from the incompressible Euler system (see e.g. or from the 

incompressible Navier-Stokes system (see for instance pi71[TTl[T71l24j ) by several 
techniques (e.g. by direct integration or asymptotic analysis). We adapt here 
the derivation in [31 3] to get a new unidirectional shallow water model. We start 
from the 3D incompressible Euler equations where we neglect the acceleration 
following the y-axis supposing the existence of a privileged main flow axis. We 
write then the Euler equations in the local Serret-Frenet reference frame in order 
to take into account the local effects produced by the changes of section and the 
slope variation. Then we derive a shallow water model by a formal asymptotic 
analysis (done in Subsection 12. 3p . 

2.1 Incompressible Euler equations and framework 

Let us consider the cartesian reference frame (0,i,j,k). In the corresponding 
coordinate system {x,y, z), the 3D incompressible Euler system writes: 



where V(t, x, y, z) denotes the velocity with components (m, v,w), P — p{t, x, y, z)!^ 
is the isotropic pressure tensor, po the density of the fluid at atmospheric pres- 
sure po and F is the exterior strength of gravity. 

We close classically System ((TJ using a kinematic law for the evolution of the 
free surface: any free surface particle is advected by the fluid velocity U and on 
the wet boundary, we assume the no-leak condition U.riwb = where n-wb is 
the outward unit normal vector to the wet part of the pipe (see FiG. We set 
the pressure P to at the free surface. 

We define the domain rii?(i) of the flow at time t as the union of sections 
fl{t, x) (assumed to be simply connected compact sets) orthogonal to some plane 
curve C lying in (O, i, k) to follow the privileged main flow axis. We choose the 
parametrization (x,0,b(x)) in the cartesian reference frame (0,i,j,k) where k 
follows the vertical direction; b{x) is then the elevation of the point lu{x, 0, b{x)) 
over the plane (0,i,j) (see FiG. [1]). 

We define a local reference frame as follows: we introduce the curvilinear variable 
defined by: 



where xo is an arbitrary abscissa. We set Y = y and we denote by Z the B- 
coordinate of any fluid particle M in the Serret-Frenet reference frame (T, N, B) 



div(/?o U) = 
9t(poU)+poU- V(poU)+VP - F 



(1) 
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at point ui{x,0,b{x)) with T the tangent vector N, the normal and B the bi- 
normal vector (see FiG. [T]and FiG. [3] for the notations). B is normal to C in 
the vertical plane [O, i, k). 

Then, at each point ui, n{t,X) is defined by the set: 

{iy,Z)eR'';Z e[-RiX),-R{X) + H{t,X)],ye[a{X,Z),P{X,Z)]} 

where R{X) denotes the radius, H{t, X) the physical water height at section 
Q{t,X). We denote a{X,Z) (respectively P{X,Z)) Y-coordinate of the left 
(respectively right) boundary point of the domain at altitude Z, —R{X) < Z < 
R{X) (see Fig. E]). We denote also -R{X) + H{t,X) by h(t,X) which is the 
Z-coordinate of the water level. 



z = R{X) 



C : z = h{X) 




e{x) 



Figure 1: Geometric characteristics of the domain 
Mixed flow: free surface and pressurized 




Figure 2: Outward unit normal n^b 7^ n (except for uniform pipes) 
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Figure 3: Cross-section il{t, X) of the domain at point uj in the free surface case 

In the sequel, we will use a curvilinear map which will be an admissible trans- 
formation under the geometrical hypothesis on the domain: 

(H) Let TZ{x) be the algebraic curvature radius of the plane curve x i-^ {x, 0, b{x)). 
We assume that: 

Vx e np, |7^(x)| > R{x). 

2.2 Incompressible Euler model in the curvilinear coordi- 
nates 

Following the work in [3l |4] , we write System ([T]) in the Serret-Frenet reference 
frame (T, N,B) at point uj{x,0,b{x)) by the transformation T : {x,y,z) — > 
{X, Y, Z) using the divergence chain rule lemma that we recall here: 

Lemma 2.1 Let {X,Y,Z) i-^ T{X,Y, Z) = {x,y,z) be a diffeomorphism 
and 

= (x,Y.zyl' ib,e jacobian matrix of the transformation with determinant 

J. 

Then, for any vector field $, one has: 

>/div(^,j^^2)$ = div(x,Y.z){J-A^) , 
and, for any scalar function f , one has: 

where A* stands for the transpose of the matrix A. 

Let ([/, V, Wy be the components of the velocity vector in the (X, Y, Z) coor- 
dinates defined as {U, V, WY = Q{u, v, w)* where O is the matrix 

(cos 9 sin 9 
1 
— sin 9 cos 9 
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where we denote by 6{x) the angle (i,T) in the (i,k) plane. 

Using Lemma 12.11 the incompressible Euler system in the variables {X, Y, Z) 

reads: 

' dx{poU)+dY{JpoV) + dziJpoW) = 

dt{JpoU) + dx{poU^)+dYiJpoUV) + dziJpoUW)+dxP = Gi 

dtiJpoV) + dxipoUV) + dY{JpoV^) + dz{JpoVW) + dY{Jp) = o 

dtiJpoW) + dxipoUW)+dYiJpoVW)+dziJpoW^) + Jdzip) = G2 

(2) 

where J{X, Y, Z) = 1 — Z6'{X) is the determinant of the transformation and 

Gi = po UWe'iX) - Jgpo sine, G2 = -po U^9'{X) ~ Jgpo cose. 

The interested reader can find the details of the calculus in 3 . We have denoted 
by /' the derivative with respect to the space variable X of any function f{X). 
On the wet boundary, the no-leak condition reads: 

{U,V,Wy .n^b ^ . (3) 

Remark 2.1 Notice that k{X) — e'{X) is the algebraic curvature of the axis 
at point uj{X,0,b{X)) and the function J{X,Y, Z) = 1 — Zk{X) depends only 
on the variables X,Z. Moreover, under the hypothesis {H), we have J > in 
flp. Consequently, T defines a diffeomorphism and thus the performed trans- 
formation is admissible. 

2.3 Formal derivation of the FS-model for free surface 
flows 

In this section, we perform a formal asymptotic analysis on System ([2]). Ac- 
cording to the work in [31 [T71 HH , the shallow water equations can be obtained 
from the incompressible Navier-Stokes equations with particular boundary con- 
ditions. Here, we perform this analysis directly on the incompressible Euler 
system in order to get J = 1 + 0{e) for some small parameter e. 
Let us introduce the usual small parameter e = H/L where H (the height) 
and L (the length) are two characteristics dimensions along the B and T axis 
respectively. Moreover, we assume that the characteristic dimension along the 
j axis is the same as for the k axis. We introduce the other characteristics 
dimensions T, P,U ,V ,W for time, pressure and velocity respectively and the 
dimensionless quantities as follows: 

U = U/U, V = eV/U, W = eW/U, 

X = X/L, Y = Y/H, Z = Z/H, p = p/P, 9^e,p^ po. 

In the sequel, we set P = and L = TU (i.e. we only consider laminar flows). 
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Under these hypotheses, we have J(X, Y,Z) = 1 — eZ6'{X). Thus, the rescaled 
System ([2|) reads: 

dj^U + dy{JV) -{^^{JW) = 
d^{JU) + di^{U^) + dy{JUV) + d^{JUW) + d^P = d 
(d^{JV) + d^{UV) + dY{JV^) + d^{JVW)^ +dY{Jp} = (4) 
+ dj^iUW) + dy{JVW) + %(JW2)^ 



where 



G2 = -eU^p{X) 2 + «'«(^) 



Fr M = , is the Froude number alone the T axis and the B or N axis 
yfgM 

where M is any generic variable equal to L or H. 
Formally, when e vanishes, System (|4]) reduces to: 

d^ij + dy{v) + d~^{w) = (5) 



sin0 



F, 



^{cosey (6) 



Let us introduce the conservative variables A(t,X) and X) representing 
respectively the wet area and the discharge defined as: 

A{t,X)^ dYdZ, Q{t,X) ^ A{t,X)U 

Jn{t,x) 

where U is the mean value of the velocity : 

We integrate the preceding system dSHHEl) along the cross-section with the 
approximation U U and UV « U V. Then, returning to the physical 

variables, the free surface model, that we call FS-model, reads: 

dtA + dxQ = 



dtQ + dx { + gh {X, A) cos 9j = gh {X, A) cos 9 - gA sin 9 (8) 

-gA'Z{X, A) (cos 6*)' 
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where Ii {X, A) and I2 {X, A) are respectively the classical term of hydrostatic 
pressure and the pressure source term defined by: 

/h ph 
{h- Z)adZ and hiX, A) = {h-Z)dxcrdZ 
-B. J-R 

which are obtained from the integration of the pressure term dj^p in Equation 
(O with p — p{h{t,X) — Z)cos9 (obtained from equation ([7])). 
In these formulas (y{X, Z) is the width of the cross- section at position X and at 
height Z. The additional term 'Z{X,A) is defined by {h{A) - Ii{X,A)/A). It 
is the ^-coordinate of the center of mass: 



/ 

Jim 



ZdYdZ 

(t,X) 
h(t,X) i'P(X,Z) 



ZdYdZ 



-R{X) Ja(X.Z) 
h(t,X) 



Z a{X, Z) dZ 

-R(X) 



>-R(X) 

In System ([5]), we may add a friction term —p^gSfT to take into account the 
dissipation of energy. We have chosen this term 5/ as the one given by the 
Manning-Strickler law (see e.g. [55]): 

Sf{A,U) = K{A)U\U\. 

The term K{A) is defined by: K{A) = — — ^ , > Q is the Strickler 

Kg Rh [A) I 

coefficient of roughness depending on the material, Rhi^A) — A/P^ is the hy- 
draulic radius and P,„ is the perimeter of the wet surface area (length of the 
part of the channel's section in contact with the water). 



3 Formal derivation of the P-model for pressur- 
ized flows 

In this section, we present a new set of unidirectional shallow water like equa- 
tions to describe pressurized flows in closed non uniform water pipes. This model 
is constructed to be coupled in natural way with the FS-model ([8]). Starting 
from the 3D compressible Euler equations in cartesian coordinates, 

dtp + div(pU) - 0, (9) 

dt{plJ) + div{pV(E)lJ) + Wp^F, (10) 

where U(t, x, y, z) and p{t^ x, y, z)) denotes the velocity with components (u, v, w) 
and the density respectively. p(t, x, y, z) is the scalar pressure and F the exterior 
strength of gravity. 
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We define the pressurized domain of the flow as the continuous extension of 
ilp (see Subsection 12. ip defined by some plane curve C with parametrization 
{x, 0, b{x)) in the cartesian reference frame (O, i, j, k); we recall that b{x) is then 
the elevation of the point lo over the plane (0,i,j) (see FiG. [1]). The curve 
may be, for instance, the axis spanned by the center of mass of each orthogonal 
section Q{x) to the main mean fiow axis, particularly in the case of a piecewise 
cone-shaped pipe. Notice that we consider only the case of infinitely rigid pipes, 
thus the sections O = fl{x) are only a;-dependent. 

We then write Equations ([ WTU)) in the {X, Y, Z) coordinates introduced in Sub- 
section [5?T] As we want a unidirectional model, we suppose that the mean flow 
follows the X-axis. To this end, we neglect the second and third equation for 
the conservation of the momentum. 

By a straightforward computation, the mass and the first momentum conserva- 
tion equation in the {X, Y, Z) coordinates reads: 



dt{Jp) + dx{pU) + dY{pJV) + dz{pJW) = 

dtiJpU) + dxipU^) + dripJUV^) + dz{pJUW) + dxp ^^^^ 

= -pjg sin e + pUW {cos 6)' 

Applying the same asymptotic analysis developed in Subsection 12.31 Equations 
PfTUl) read: 



dt{p) + dx{pU) + dY{pV) + dz{pW) = 

dt{.pV) + dx^pV"^) + dY{pUV) -f dz{pUW) -f dxv = -p.gsin6i 

— .gZ(cos( 

We choose the linearized pressure law: 



(12) 



(see e.g. [29l[30]) in which po represents the density of the fluid at atmospheric 
pressure poj Pa is some function set to zero and (3 the water compressibility 
coefficient (equal to m? .N~^ in practice). The sonic speed is then 

given by c = l/\//9po ^md thus c w 1400 m.s~^. 

For m e 9f2, n = ■; — r is the outward unit vector at the point m in the il-plane 
|m| 

and m stands for the vector Lom (as displayed on FiG. [3]). 

Following the section-averaging method performed in Subsection 12.31 we inte- 
grate System (IT^ over the cross-section $7. Noting the averaged values over J7 by 
the overlined letters (except Z), and using the approximations pU ~ p?7, pU'^ ~ 
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2 

pU the shallow water like equations read: 

dtipS) + dx(.pSU) = [ piUdxni-Y).nds (14) 

Jon 

dt ipSU) + dx {pSU + c^-pS) = -gpS sin 9 + c^pS" 

- gpSZ{cosey (15) 

+ / pU {Udxm-V) .nds 
Jan 

where V — {V,WY is the velocity in the (N,B)-plane. We denote by S the 
area of the cross-section of the pipe at position X. 

The integral terms appearing in (jl4p and IjlSp vanish, as the pipe is infinitely 
rigid, i.e. fl = il{X) (see [8^ for the dilatable case). It follows the non- 
penetration condition (see FiG. SJ: 



•nwb = . 



Omitting the overlined letters (except Z), we introduce the conservative vari- 
ables 

A ~ —S the FS equivalent wet area (16) 
Po 

Q — AU the FS equivalent discharge . (17) 
and dividing Equations (fT4| -(|15 p by po we get: 
dt{A)+dx{Q) = 

dt{Q)+dx (^ + c^a] = -gAsinO- gAZ{X,S){cosey (18) 




A 



-'4 



As introduced previously for the FS-model in Section (|2.3p . we may introduce 
the friction term —pgSf T given by the Manning-Strickler law (see e.g. [^1: 

Sf{S,U)=K{S)U\U\ 

where K{S) is defined by: K{S) = ^(Q\i/^ ' i^s > is the Strickler co- 

efiicient of roughness depending on the material and Rh{S) = S/Pm is the 
hydraulic radius where Pm is the perimeter of the wet surface area (length of 
the part of the channel's section in contact with the water, equal to 2 tt i? in the 
case of circular pipe). 
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This choice of variables is motivated by the fact that this system is for- 
mally close to the FS-model (|8]) where the terms gli (X, A) cos 0, g/2 {X, A) cos 0, 
— S' — 

Z{X, A) are respectively the counterparts of c^A, c^A—, Z{X, S) in System (US]). 

o 

Let us remark that the term Z is continuous through the change of state (pres- 
surized to free surface or free surface to pressurized state) when the same curve 
plane is chosen (in practice, the main axis of the pipe). Then, we are moti- 
vated to connect "continuously" System ([8]) and (fT8|) through transition points 
(through the change of state) by defining a continuous pressure law. It leads to 
a "natural" coupling between the pressurized and free surface model as we will 
see in Section |4l 

4 The PFS-model 

The formulations of the FS-model and P-model ([T5)) are very close to each 
other. The main difference comes from the pressure law. In order to build a 
coupling between the two models, we have to define a pressure that ensures its 
continuity through transition points in the same spirit of [7] ■ As pointed out in 
the previous section, we will use the common couple of unknowns {A, Q) and 
the same plane curve C (see Remark 14. ip to get a continuous model for mixed 
flows. 

Remark 4.1 The plane curve with parametrization (x, 0, b{x)) is chosen as the 
main pipe axis in the axisymmetric case. Actually this choice is the more conve- 
nient for pressurized flows while the bottom line is adapted to free surface flows. 
Thus we must assume small variations of the section (S' small) or equivalently 
small angle ip as displayed on FiG. [7} 




Figure 4: Some restriction concerning the geometric domain. 
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We introduce a state indicator E (see FiG. H} such that: 

^ _ f 1 if the state is pressurized: (p ^ po) ,^g^ 
[0 if the state is free surface: {p — po) ' 

Next, we define the physical wet area S by: 

s = s(Ai.).{^ ;[ (20) 

and a modified pressure law (see FiG. H]) which ensures its continuity through 
the change of state by: 

p{X, A,E) ^c^{A~S)+ gh {X, S) cos 6*. (21) 

Remark 4.2 

• Indeed, when a change of state occurs we have: 

hm p{X, A, E) = hni p{X, A, E) = ghlX, S) cos(e) 

A<S A>S 

which ensures the continuity of the pressure. 

• The flux gradient F is discontinuous through the change of state since 

Qp Q gp 

— {A, Q, 0) = 9qjIi{X, A) cose ^ c^ = —{A, Q, 1). 

FinaUy, from the P-model (HH]), the FS-model ©, the definition of E (HH), the 
definition of S d^O]) and the pressure law (pij) . the PFS-modcl for unsteady 
mixed flows can be simply expressed into a single formulation as: 

dt{A) + dx{Q) = 

dt{Q)+dx [^+p{X,A,E)] = -gAb' + Pr{X,A,E) 

^ J (22) 

-G{X,A,E) ^ ' 

-K{X,A,E)9^ 

where Pr, and G denotes respectively the friction, the pressure source and 
the geometry source term defined as follows: 

Pr{X,A,E) = c^ (^^-l^ S' + gl2{X,S)cose 

f-H(S 
l-R{X) 



with/2(X,S)= / {n{S) - Z)dx<7{X,Z)dZ, 

J-R(X) 



G{X,A,E) = gAZ{X,S){cosey, 
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and b' stands for sin9{X). H represents the Z-coordinate of the water level: 

Remark 4.3 (Both models are recovered) Setting S{A, E) — A in System 
ll2Si) . we obtain obviously the free surface model 0). For all pressurized states, 
when S{A, E) = 5*, the pressure law i21\) reads, for instance, in the case of 
circular pipe: 

c^{A~ S) + gIiiX,S) cose = c^{A - S)+gTrR^ cosB 

which is not exactly the pressure law of the P-model il8\) . Indeed, the derivation 
of the P-model is done with the linearized pressure law il3\) (see Section\3^ with 
Pa — 0. Thus, the property of the continuity of models il8\) -l[^} through a change 
of state is obtained if and only if Pa is chosen as gIi{X, S) cosO which is the 
hydrostatic pressure corresponding to a full section. 




Sl,{X.A) cose 



S = A 




Piezometric level over the bottom 



Figure 5: Free surface state A, 0) — g Ii{X, A) cos9 (top), pressurized state 
with overpressure p{x, A,l) > (bottom left), pressurized state with depression 
p{x,A, 1) < (bottom right). 
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The PFS-model satisfies the following properties: 
Theorem 4.1 

1. The right eigenvalues of System \2^) are given by: 

\- = U - c{A, E), X+ = U + ciA, E) 



with c{A, E)^ ^ ]j3 T{A) ^ ^ ^ , where T{A) is the width 

c if E^l 

of the free surface (see FlG. 

Then, System is strictly hyperbolic on the set: 



{A{t,X)>0} . 

2. For smooth solutions, the mean velocity U — Q/A satisfies 

= -gK{X,A,E)U\U\ 

[/2 



dtU + dx\^^+c^ HA/S) + gn{S) cos + gb j ^^4) 



The quantity — — h ln(A/ S) + gTl{S) cos 9 + gb is called the total head. 
3. The still water steady state reads: 

u = and ln(yl/5') + gn{S) cos O + gb^Q. (25) 

4-. It admits a mathematical entropy 

n2 _ 

£{A,Q,E) = + Ahi[A/ S) + S + gAZ{X, S) co^e + gAb (26) 

which satisfies the entropy relation for smooth solutions 

dt£ + dx(^{£ +p{X,A,E)) u'j = -gAK{X, A, E)U^\U\ ^0. (27) 

Notice that the total head and £ are defined continuously through the transition 
points. 

Remark 4.4 The term AZ{X, A){cos9y is also called "corrective term" since 
it allows to write the Equations ^24^ and {21^ with i26\) . 

Proof of Theorem I4.lt the results (l24l) and l(27|) are obtained in a classi- 
cal way. Indeed, Equation ([M)) is obtained by subtracting the result of the 
multiplication of the mass equation by U to the momentum equation. Then 
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multiplying the mass equation by f — + ln(A/S) + g'H{S) cos 6* + gb \ and 
adding the result of the multiplication of Equation (IMl) by Q, we get: 

dt + c^A ln(A/S) + c^S + gAZ{X, S) cos 9 + gAh 
+dx ( + ln(A/S) +c^S + gAZ{X, S) cos 6* + gAfo + p{X, A, E)J U 
- l) = -9AK{X, A, E)U^\U\ . 

We see that the term ^— — 1^ dtS is identically since we have S = A 

when the flow is free surface whereas S = S{X) when the flow is pressurized. 
Moreover, from the last inequality, when S = A, we have the classical entropy 
inequahty (see [Illl]) with £: 

£{A, Q, £:) = ^ + gAZ{X, A) cos + gAb 
while in the pressurized case, it is: 

£{A, Q, £;) = ^ + c^A \n{A/S) + c^S + gAb. 
Finally, the entropy for the PFS-model reads: 

E{A, Q,E) = -^ + c^A \n{A/S) + c^S + gAZ{X, S) cos 9 + gAb. 

Let us remark that the term c^S makes £ continuous through transition points 
and it permits also to write the entropy flux under the classical form (£ +p)U. 



5 Finite volume discretisation 

In this section, we adapt the VFRoe scheme described in [TH [111 [7] . The new 
terms appearing in the PFS-model related to the curvature and the section vari- 
ation are upwinded in the same spirit of [7] . The numerical scheme is adapted to 
discontinuities of the flux gradient occurring in the treatment of the transitions 
between free surface and pressurized states. 

5.1 Discretisation of the space domain 

The spatial domain is a pipe of length L. The main axis of the pipe is divided 
in cells nii = [Xi_i/2, X^^l], 1 < i < N . At" denotes the timestep and we set 

tn+l = tn + At". 
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The discrete unknowns are UJ' 



For the sake of simpHcity, the 



boundary conditions are not treated (the interested reader can find this treat- 
ment in details in [J). 



5.2 Explicit first order VFRoe scheme 

We propose to extend the finite volume discretisation [7] to the PFS-model using 
the upwinding of the new source terms: the curvature and section variation of 
the pipe. In what follows, we do not write the E dependency. 
First, following Leroux et al. [111 [55] we use piecewise constant functions to 
approximate b {b'{X) = s'md{X)) as well as the term cos 61 and the cross section 
area S. Adding the equations dtZ — 0, dtC0s9 — and dtS — 0, the PFS- 
model can be written under a non-conservative form with the variable W = 
{b,cose,S,A,QY: 



dtW + dxF{X, W) + B{X, W)dx'W = TS(W) 



(28) 



where 



F(X, W) 



, TS{W) = 



-p{X,A) J 



\ -9K{X,S) 



Q\Q\ 
A 



and 



B{X,W) = 



/ 





0_ 

V gA gAZ 



\ 









-c2(A/S-l)-X(X,W) 0/ 



where we have written the pressure source term due to the geometry gl2 {X, S) cos(0) 
as X(X, W)S". For instance, for a circular cross-section pipe we have: 



„(S).,.csi„(||>) 



Let WT* be an approximation of the mean value of W on the mesh rrii at time 
tn- Since the values of &, cos 0,5* are known, integrating the above equations 
over ]Xi_i/2i X^j^i [x [i„, [, we can write a Finite Volume scheme as follows: 



= Wr-a,(F(W*+i/2(0~,Wr,WIVi))-F(W*_i/,(0+,W^'_i,Wr)) 



+T5(W7) 



(29) 
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with ai = . 

hi 

Wi+i/2(f — Wi, Wi+i) is the exact or an approximate solution to the 
Riemann problem at interface Xi^i/2 associated to the left and right states 
Wi and Wi+i. Let us also remark that the term B{X,W) does not appear 
explicitly in this formulation since b', (cosO)' and S' are null on ]Xj_i/27 -^j+i [ 
but contributes to the computation of the numerical flux. 

The computation of the interface quantities W*_|_j^/2(0*i Wj+i) will depend 
on two types of interfaces located at the point X^^i : the first one is a non 
transition point, when the flow on both sides of the interface is of the same 
type. The second one is a transition point, when the flow changes of type 
through the interface. We recall the approach used in [7 and adapt it here to 
the new terms. According to the type of interface, we have to solve two different 
linearised Riemann problems. 



5.2.1 The Case of a non transition point 

Expanding the term dxF{X,W) in the non-conservative equations ((28)) . the 
quasilinear formulation of the PFS-modcl (|22|) reads: 

(30) 



dtW + D{W) dxW = T5(W) 

with D the convection matrix defined by 

/ 







\ gA gAn{S) l'(W) c'^(W) ~ 

where *(W) = gSdsHiS) coa0 - c^(W)- and u = Q/A denotes the speed of 



D(W) 



\ 




1 

2u J 



(31) 



the water. c(W) is then c for the pressurized flow or 
surface flow. 



A 



TiA) 



cos 9 for the free 



Remark 5.1 Let us remark that, since dxhiX, A) ~ l2{X, A) + dAli{A)dxA, 
the pressure source term I2 does not appear in Equation fi30(l . 

To compute the interface quantities denoted by [AM, QM) for the left hand 
side and [AP, QP) for the right hand side (see Figure [6] below) , we solve the 
following linearised Riemann problem: 



W 



D dxW = 



W; = ibi,cosei,Si,Ai,Qiy 

Wr = {br,C0s6r,Sr,Ar,QrY 



if 
if 



X <0 
x>0 



(32) 

with (W/, W^) = (Wi, W,+i) and D = D(W/, W^) = D{W) where W is some 
approximate state of the left W; and the right state. 
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Remark 5.2 We will see in Section\^ that the classical approximation D{ W) 

of the Roe matrix DRoei Wi, Wr) = / D{Wr + (1 - s){Wi - Wr)) ds defined 

Jo 



by D = D{W) = D 



Wl + Wr 



is not a suitable choice to preserve the still 

water steady state. However, we propose in Section\E\ a new approximation of 
D which maintains it perfectly. 

We have then W^*(0+, W,, W^)^= {br, cosOr, Sr, AP, QPf . _ 
The eigenvalues of the matrix D are Ai = 0, A2 = 0, A3 = 0, A4 = u — c(W), 
A5 = u + c(W) and the associated right eigenvectors: 



/ C2(W) 





ri(W) 



0^ 

-gA 




r2(W) 



/ *(W) \ 

0^ 

-gA 





r4{W) 



r,(W) = 



\u- c(W) / 



r3(W) 



-1 






V u + c{W) J 



(1),(2),(3) (1),(2),(3) (1),(2),(3) 




W| Wr Wl Wl Wr 

U>C U<-C -C<U<C 



Figure 6: Solution of the Riemann problem (|32|) . The lines (i) corresponds to 
the characteristic lines X/t = Ai, for i = 1, . . . , 5 . 

We denote P the transition matrix associated to the right eigenvectors of D and 
P~^ its inverse. Setting [[W]] = W^ — Wj, the solution of the Riemann problem 
are constant states connected by shocks propagating along the characteristic 
lines X/t — A^. The jump associated to the eigenvectors is then equal to 
{P^^ ri. In particular, the discharge is continuous through the line X/t = 

since the fifth component of vectors ri, 7-2 and r^ are null. It follows that the 
equation associated to the wet area A remains conservative. 
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Thus, for instance in the subcritical case (when — c(W) <u < c(W)), we have: 

AM = Ai + ^-^4^ ^[ + }L±^!^ - Ai) \^ {Qr - Qi) 

2c(W)(c(W)-S) 2c(W) 2c(W) 

QM = QP = Q^-^i,l + ^—^{A,-A,)-^-^^(Q.-Q^) 
2c(W) 2c(W) 2c(W) 

~2-.(W)2 



— c 



where V'f is the upwinded source term br — bi+'H{^){cos6r — cos6i) + ^{yf){Sr — 
Si). 

Remark 5.3 The friction term can also be upwinded in the same way. Writing 
the friction term under a conservative form 

8. r Kis,S)^^m¥^ds 

Jxo 



(for some arbitrary Xq) allows us to write the "static" slope b as a "dynamic" 
one as follows: 

Jx A^t,s) 

that we denote again b. Thus, the upwinding of the dynamic slope bi+i — bi 
reads: 



1 j_Q\Q\ \ r^'+^ 1 [ Q\Q\ \ 



which is equal to: 

since A and Q are constant on each cells. 

The terminology "dynamic" and "static" slope is used since one is {t, x)-dependent 
while the other is only x-dependent. 



5.2.2 Case of transition point 

In the case of a transition point, we assume that the propagation of the interface 
(pressurized-free surface or free surface-pressurized) has a constant speed w 
during a time step. The half hne x = is then the discontinuity line of 

D{Wl,Wr). 

Let us now consider U = {A~ ,Q~) and U'*' = (5+) the (unknown) states 
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respectively on the left and on the right hand side of the line x = wt with 

Q+ - Q- _ , , , 

w — — — -— . Both states U; and U (resp. U,- and U ) correspond to 

the same type of flow. Thus it makes sense to define the averaged matrices in 
each zone as follows: 

• ioi X < wt, we set Di — Z3(W/,W,.) = Z3(W/) for some approximation 
W/ which connects the state W; and W~ (see Remark l5.2p . 

• for X > wt, we set Dr ~ Z)(W;,W.r) D{Wr) for some approximation 
W; which connects the state and (see Remark l5.2p . 

Then we formally solve two Riemann problems and use the Rankine-Hugoniot 
jump conditions through the line x = wt which writes: 

Q+-Q- = w{A+-A-) (33) 
F5{A+,Q+)-F5{A-,Q-) = w{Q+-Q-) (34) 

with F5{A,Q) = +p{X,A). According to the left (U", UM) and right 

unknowns (U"*", UP ) at the interface Xi+1/2 and the sign of the speed w, we 
have to deal with four cases: 

• pressure state propagating downstream, 

• pressure state propagating upstream, 

• free surface state propagating downstream, 

• free surface state propagating upstream. 

We can next consider two couples of "twin cases" : pressure state is propagat- 
ing downstream (or upstream) as shown in the figure [7] and free surface state 
propagating downstream (or upstream) as shown in the figure [51 Moreover, for 
all existing transition case, the upwindcd altitude term br — 6; in are replaced 
by V-r- 

Pressure state propagating downstream: it is the case when on the left 
hand side of the line ^ — wt, we have a pressurized flow and on the right hand 
side we have a free surface flow: the speed w of the transition point being 
positive. Following Song [2H] (see also [H]), an equivalent stationary hydraulic 
jump must occur from a supercritical to a subcritical condition and thus the 
characteristics speed satisfies the inequalities: 

Ur + c(W)r < W < Ul + C 

where c is the sound speed for the pressure flow, u;, Ur, and c(W)r are defined 
by the same formula obtained in the case of a non transition point but according 
to Di and Dr- 
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Figure 7: Pressure state propagating downstream. 



Therefore, only the characteristic Hues drawn with solid lines are taken into 
account. Indeed they are related to incoming waves with respect to the cor- 
responding space-time area — oo < ^ < w. Conversely, the dotted line ^ = 
Ur — c(W)r, for instance, related to the free surface zone but drawn in the area 
of pressurized flow is a "ghost wave" and is not considered. Thus = and 
U;, U~ are connected through the jumps across the characteristics ^ = and 
£^ — ui — c. Eliminating w in the Rankine-Hugoniot jump relations (|33|) -p4 p . 
we get U~ as the solution to the nonlinear system: 



{F5{Ar,Qr)-F5{A-,Q-)) = 



-Qi-{A -Ai){ui-c) + 



Finally, we obtain 



' AP 


= A- 


QM 


= Q- 


< QP 


= Q- 


AM 


= AP 



{Ar-A-) 
C + Ul 



uf — c^' 



= 



(35) 
(36) 



Free surface state propagating downstream: on the left hand side of the 
line £^ — wt we have a free surface flow while on the right hand side, we have a 
pressurized flow (the speed w of the transition point being positive) . Following 
Song again, the characteristic speed satisfies the inequalities: 

Ul + c(W)i < W <Ur + C 
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Figure 8: Free surface state propagating downstream 



There are two incoming characteristic hues with respect to the free surface area 
— oo < ^ < w (actually three with ^ = 0) and they can connect the given left 
state U; with any arbitrary free surface state UM. Thus only one characteristic 
line — Ur+c) gives any information (it is the equation ([37]) above) as an incom- 
ing characteristic line with respect to the pressurized zone w < £, < +oo. From 
the jump relations through the characteristic ^ = 0, and after the elimination 
of w in the Rankine-Hugoniot jump relations (|33p . (|34p we get another equation, 
namely Equation ([55)) above. It remains to close the system of four unknowns 
{A~ , Q~ , , Q"*")- Firstly, we use a jump relation across the transition point 

(with speed w) for the total head ^E* = '"'^^^^ ( "5" ) + g'HiA) cos9 + gb^ from 



Equation ()24p . which writes: 

w (ii+ - u~) . 

Finally, we use the relation: 

Qr — Qi 



W = Wpred with Wpred 



Ar - Ai 
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We have then to solve the nonlinear system: 

{Qr - 0+) = {Ar - A+) {Ur + c) (37) 

(Q+ - Q~) {Qr - Qi) = [Ar - Ai) (F2(A+, 0+) - F2(A", Q-)) (38) 
' In {A+) + g cos 6 n{A+) - -^^-^ - In {A-) ^ g cos 9 n{A-) 



- [a^ - A^J 

{Qr - Qi) {A+ - A-) = (Q+ - Q-) (A, - Ai) (40) 

The states UM et UP are then obtained by the following identities: 

gAiiiT ui +c(W)i , ^ , , 1 , 

AM = Ai + 'Jl' + {A- - Ai) ^ {Q- - Qi) 

2c{W)i{c{W)i^ui) 2c{W)i 2c{W)i 

AP = AM' 



«?-c(W)f 

gAiVi 



QM = QP = QMP = Qi + 



2c(W)j 



_^ iIf-c(W)^ (^-_^^)„ ^.-c(W). 
2c(W)z 2c(W), 

Finally, the update state A^'^^ , Q"^^ are obtained by the same relation as in 
the case of a non transition point. 

Using equations we update the values of A"^^, Q"^^ with a standard sta- 
bility condition of Courant-Fricdrich-Lcvy controlling the time step size Ai" . 

5.2.3 Updating the state of the flow E va a cell. 

To update the state E in the cell (see FiG. [9]), we use a discrete version 
of the state indicator E defined by equal to 1 for a pressurized flow and 
otherwise. Following [7], after the computation of the wet area ^""""^ we predict 
the state of the flow in the cell rrii by the following criterion: 

• if Ef = then : 

if A^+i < S^ then = 0, else E'l^^ = 1, 

• \iEf ^l: 

if A1+^ > S, then E^+^ = 1, else = £;f_i • E^f+j. 

Indeed, if A'^'^^ > Si it is clear that the state of the flow in the cell nii becomes 
pressurized, on the other hand if A^^^ < Si in a, mesh previously pressurized. 
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we do not know a priori if the new state is free surface (p = po and the value of 
the wetted area is less than Si) or pressurized (in depression, with p < po and 
the value of the wetted area is equal to Sf. see Remark and FiG. [TU)) . 
So far as we do not take into account complex phenomena such that entrapment 
of air pockets or cavitation and keeping in mind that the CFL condition ensures 
that a transition point crosses at most one mesh at each time step, we postulate 
that: 

1. if the state of the flow in the cell rrii is free surface at time tn, its state at 
time tn+i is only determined by the value of A^^^ and it cannot become 
in depression. 

2. if the state of the flow in the cell is pressurized at time i„ and if 
A""*"^ < Si, it becomes free surface if and only if at least one adjacent 
cell was free surface at time t„. This is exactly the discrete version of the 

continuous — criterion in Remark 15.41 and displayed on FiG. [TUl 



yes, if Aj < A 




Figure 9: Update of the state E"~^^ of the mesh m^. 



Remark 5.4 As we do not take into account complex phenomena such that 
entrapment of air pockets, each connected component of the pressurized area 
is simply connected (see FlG. \1U\) . Moreover, for each depression area D, its 
closure D is a strict subset of the pressurized set. It follows that when A < S 
on each pressurized area, we observe a depression as displayed on FlG. \10l 
Moreover, we may also use a visual depression indicator given by the function 

— .• the case S ^ A corresponds to a free surface state while S > S to an 
overpressure state and S < S to a depression state. On FlG. \1(A we draw the 
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behavior of the interface 

— at fixed time t^ . 
o 



w in the (X, t) -plane and the graph of the function 



t = to 



t = h 



t = tn 



t = t3 



t = ti 



t = t4 
t = t3 
t = t2 
t = tl 
t = to 



s 





dX 

w = —r- 
dt 




Dressu 


rise 




e 


dX 

w — —— 
dt 








rest 






















free surfaa 


i zon 


e^~^ 






free 


surface zone 
















1 FS state: S = A 


overpressure 
^ > S 
S = S 

p(X,A,l)>0 




overpressure 
A > S 

s = s 

p(X,A, 1)>0 


X 

FS state: S = A 













depression 
^ < S 

s = s 

p{X,A,l)<Q 



X 



Figure 10: — depression indicator. 
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6 Remarks on still water steady state: an ex- 
actly well balanced scheme 

This section is devoted to the construction of an exactly well-balanced scheme in 
the sense that it maintains perfectly the still water steady states. This scheme, 
noted EWBS, is obtained by a suitable definition of the convection matrix. 
The numerical approximation of the PFS-model (l?T|) reads: 

At 



= ^r-ii^(Q^Vi/2-Qr-i/2) (41) 



Qr' = Qr-^(^5(AMf+i/2,Qr+i/2)-^^5(Ai^ii/2,Qr-i/2)) (42) 

g2 

where Qi±i/2 stands for QMP,±i/2 andF5(A, Q) ^ -^ + c^{A-S)+g Ii{X,S). 
For instance, in the subcritical case, the interface quantities reads: 

a A" 

U", , /„ + Ci + l/o 1 



Qr+1/2 = Q^-L^^|+i+ '+y^ '+^/^ (A^Vi-^r) (43) 



2 

^-^1+1/2 — ^^'"'1+1/2 + ~„,2 ~2 

"1+1/2 ^i+1/2 

where the upwinded source term reads: 

6,+i - 6, +7^(Sr+i/2)(cos0,+i - COS0,) + *(Wz + l/2")(5,+i - S,) 

and 0^+1/2 stands for c\y^ ^j^-^j^) with 

Wj_,_l/2 = (^1+1/2, COs6'i_|_i/2, 5*1+1/2, ^"+1/2j '3"+1/2) 

given by 

~ _ fa; + ' _ COS Qi + cos 6*,+! 2; _ + Si+i ~„ _ Qr + <3r+i 

— ^ 1 cosy - ^ , 5 — ^ , (^, + 1/2 - ^ , 

(44) 

and the approximation of ^"_^i/2 specified. 

Starting from a discrete state (A", Q") at time t„ such that: 
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(H) let n such that: Vi, Q" = and satisfy the discrete still water steady 
state equation (according to Equation ([24]) ): 



c 



4" \ /A 



i+l 



In 1 ^ )+5H(Sr+i)cos0,+i+.95,+i = c^ln ( ^ ) +5^(8^) cos , 

(45) 



we will say that: 
Definition 6.1 

1. The numerical scheme 114^4^ - i44^ for some approximations of the terms 
^r±i/2 {^A^kq) well-balanced (also denoted by [ItA, kQ)-WB) if: 

Vi, |v4^+i-^f| = 0((max,gz/jO'='*) and |g^'+i - | = 0((max,ez /lO'^Q) 

with ItA > ^, kg > 1 the well-balanced order of the numerical scheme 
and ^4^ respectively. 

2. The numerical scheme 141M^ -144^ for some approximations of the terms 
^r±i/2 exactly well-balanced (also denoted by EWB) if: 

We will denote by {Ua, fcQ)-WBS the (fc^, kq) well-balanced scheme and EWBS 
the exactly well-balanced scheme. 

(SF) In the rest of this paper, we assume hi = AX constant, the radius R 
and b are given afhne functions, the angle 9 is constant which implies 
that the jumps across the interface Xi+1/2: — Ri+i — Ri = Ai?, 

A5i_|-i/2 — fej+i ^ bi — Ab are constant and A cos 6^+1/2 is null. 

Then under this simplified framework, wc show that: 

Theorem 6.1 

1. • The numerical scheme with the classical choice 

An , An 

^r+1/2 - (46) 

and non constant section S is not well-balanced in the sense of Defi- 
nition [KT\ (we have kg = 1). 

• For constant section and Z — 0, the numerical scheme |^JH^^[ j- |^^[ j 
with is EWB. 

2. Under a suitable choice of A^^-^^^^, the numerical scheme |^JH^^ - |^^[ j is 
EWB. 

The following section deals with the numerical scheme with A^^^^^ defined as 
([46|l where we show the first point of Theorem l6.ll The second point is studied 



in Section \6l2\ where a convenient definition of A"^^^^ leads to an EWBS. 
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6.1 Still water steady state and the classical approxima- 
tion 

The simpler choice of definition for the convection matrix D(W) is the one ob- 
tained by the approximation of the mean value of the Roe matrix Dj^oe (W; , Wj.) = 

/ D(Wr + (1 — s)(W( — Wr))ds. This approximation is given hy D = 
Jo 

I?(W) = D I ^ I that we call "the classical approximation". Thus, 



An 1 An 

Jn _ 
^j+1/2 — 9 



^ ^ 2 ^ 
defining A^^^^^ as follows: 



provides the classical approximation. But, it is not suitable to preserve the still 
water steady state: we will see that the numerical scheme (|41j|42p - (|44p with (|46p 
defines a non well-balanced scheme in the sense of Definition (|6.ip since kg = 1. 
To this end, let us assume (H) and (SF) at time i„, Equations ((43|) read for 
every i: 

/^ J" A 4" 
Vi+l/2 ~ ~"o^^ '^^i - C,+ i/2 • 

(7 A" 

S+1/2 

Denoting Qi+i/2 - (9i-i/2 by AQ^+i/z, we have: 

^Ci_l/2 Cj+1/2 

-(iiVi/2 -^1/2) 2^+1/2 
+ (H.+1/2 - c,_i/2) ir+1/2 v^+'l 



^^"+1/2 ,~ ~ X 

7y l'^i-1/2 ~ Cj_,_i/2j 
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where 



An , _ An 



i+1 



= < 



+1/2 



-Ab if = 



c2 AS-, 



2 

+ 1/2 



'i-1/2 'S'i+i/2 



I X [as., 



1 



^ CAX if i?j = ( for some constant C) 
= if iJj = 1 ( since Cj+1/2 = Ci_i/2 = c) 



|Ci+l/2 - Ci-1/2| 

Denoting tlicn 

M = max (mp (q+i/2 I?_i/2) , max (q+1/2 1^*+^ |) , max (1^+1/2 l) > c) 
and observing that 

Vi, 0{AR) = 0(A5i+i/2) = 0(AA^Vi/2) = ^(AX) 



and 

we deduce: 
It foUows that: 



0(Ci+i/2 Ci_i/2) = 0(5,+ i/2 S'j_i/2) = 0(1) , 

\AQ1^,/^\^MAx\ 



n+1 4 n 



Then, we denote i^sCAM;;,/^, Q^^^/J - ^^-1/2) by Ai^ = Ti + T2 

with 



Ti = 



and 



T2 



AM^y^ - ^^r-1/2) (Qr+1/2)' + ^AflVl/2 ((Qr+1/2)' - (Qr-1/2) 

^^r+1/2^^^-1/2 



g COS e (Ji (X,+i/2, , ^Mf+1/2) - /i (Xi+1/2, , ^f^-1/2)) if = 

c'(^^r+l/2-^^-l/2)+c'A5i+i/2 if E, = l 
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/ s 2 

where | (^Q"+i/2 



2+1 / 2 i—l / 2 




(48) 



As the term — ^ 1 is at least of order AX since 

\ ''i+l/2 '^i-l/2 / 



we have: 

It follows that : 



\AMJX^/2-APr_i/2\^OiAX). 



\Q7^' - Q7\ = OiAX) . 

Remark 6.1 For constant section S and Z ^ 0, it is easy to see that the 
numerical scheme |^ with ^JS^ is EWB. 

Although the scheme ()4m42p -(|44 |) with (|46|) have an order kq = 1 for non 
constant section, the still water steady state for the pressurized case is very well 
maintained for great value of the sonic speed c. But it is not the case for the 
free surface numerical scheme (see FiG. [T6]with AX = 10""^, 6 = 10~^, uniform 
pipe with diameter 1). We plot in FiG. [14] and FiG. [15] a still water steady state 
for two values of c for a given AX = 1 ( with h = —0.9 ): the absolute error 
obtained is 10~^ for c = 30 and 10~^ for c = 200. Indeed, let us consider, for 
the sake of simplicity, the case S = cte. At the continuous level, the still water 
steady state equation reads: 

In (A) + gRcos + gb = cte. 
With the hypothesis (H), in particular using Equation ((45|) . we write: 

A«+i=Arexp(-^A6) . 
Given AX , for great value of c, we can approximate A^^i by: 

A" 1 » A" f 1 - 4a6 



Then, replacing the right hand-side of A2_^^ in ([^5]) . we have 

\AMl]_,/,~APr_,/,\^g^:\^^o(^^^ . 

As we can see kn is always equal to 1 but the constant ^ plays the role of 

a smoothing term which helps the scheme (|4m42l) - p4)) with (l46l) to stabilize 

rapidly towards the equilibrium. Physically, c is approximatively 1400 (for a 

1 

pressurized flow without air), thus ^ w 1.9 10 . Since c(A) ^ c, this feature 
is not observed for the free surface numerical scheme. 
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6.2 An exactly well-balanced scheme 



This section is devoted to the construction of an EWBS. We have seen in Subsec- 
tion [^3] that the classical approximation of with SSI is not appropriate 
to preserve the still water steady state. Thus, we have to find a suitable defi- 
nition for A 



i±l/2 



to obtain an exactly well-balanced scheme. For this purpose, 
let us assume (SF) and start with : 

at the discrete level, the still water steady state is perfectly maintained (see 
Fig. [TT|) : it exists n such that for every i, if Q" = and Vi, 



Al: 



' A'- 



5. 



i+l 



AT. 

-<?H(Sr+i)cos0 + 56.+i =c^ln( ^ ) + gHiS^) cosO + gh 



A2: = AP-_,^„ 

A3: = Q"-i/2j 

then, for all / > n the conditions Al, A2 and A3 holds. 



Z A 



,.0-1/2 



AM, 



+ 1/2 




^io+1/2 



X 



Figure 11: Discrete representation of the mixed free surface-pressurized still 
water steady state at time t". 

The condition A2 is satisfied if and only if 
AMr,,,2 - APr^,,2 . A^r.v. + f + = . (49) 

^ \ %+l/2 ^-1-1/2 / 

The condition A3 is satisfied if and only if 

9 / ^1-1/2 ^r+1/2 ^1^' 1 ^^r+1/2 ~ . 

^Q^+l/2 = ^ \ ~ ~ ) + 7, (c,_i/2 - Ci+1/2) - . 

(50) 
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The condition Al is satisfied for pressurized flows if and only if 

Ar+i=Ar^exp(-^(A6 + Ai?cos0)) . (51) 

The condition Al is satisfied for free surface flows if and only if 

_ Kcos9-Ab 

For circular cross-section pipe, A'2_^i is computed by: 

^r+i = ^(c..+i-sin(c..+i)) (53) 

with u)i+i = 2 [ TT — arccos ( ] is the angle displayed on FiG. 



R 



Figure 12: angle lo 

Thus, the discrete still water steady state is perfectly maintained if and only if, 
(A"_j^^27 ^"+1/2) is solution of the non-linear system: 

- AA" +^( ^+^/^^^^' + ^ti/2^l-i \ 

^ \ ''1+1/2 S-1/2 / 

g j ^7-1/2 ^i+1/2 i^l^^ 1 , ^^"+1/2 ~ X 

= i } + ^ (Ci-1/2 - C,;+i/2) 

(54) 

where we have replaced the expression of Af_^i in (|49ll50p by ((5T|) for pressurized 
by ([5^ for free surface flows: 

4" (^^exp(-J(A6 + Ai?cos0)) -1^ if ^, = 



cos 9 

with T : hi-^ ^(h) = A. For circular pipe, ^ is given by 
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Finally, the numerical scheme (|41j|42p - (|44p with ^"j_]^y2 solution of the 

non linear system ((54|) defines an exactly well-balanced scheme. 

For uniform pipe and pressurized flow, the previous system simply writes: 



r+l/2 + + ^"-1/2) - 

An ^ An 

^1+1/2 — ^1-1/2 

The solution is easily obtained by: 

T„ _ ,2 4"(exp(-^A5)-l^ 

^^+1/2- 7 A5 " "7 A5 ■ ^^^^ 

Let us also remark that using the relation A" = exp A6^ , we have: 

_ c^A^m/._ c^ ^m(l-exp(^A6)) 
^^+1/2- 7 A6 " 7 A5 • ^^^^ 

It follows that A-Yi/2 expressed as the mean value of ([55| and ((56)) as 

follows: 



A-n. 

^+1/2- 5A6 



For small AX, we have 



,2 r Ar+i (1 - exp (|A6)) + 4" (exp (-^Ab) - l) 



/in I An 
An ^ ^ 



i+1/2 ~ 2 

It follows that the scheme (|4m42l) - PH) with is the zero order approximation 
of the solution given by the EWBS. 

The same analysis shows that the free surface numerical scheme with pS)) is 
also the zero order approximation of the solution given by the EWBS. 



On Fig. [TH Fig. [13 Fig. [TH we display the stiU water steady pressurized 
and free surface state computed by the EWBS and the scheme with the approx- 
imation (Uni)- The steady state with the EWBS is numerically well preserved 
while as pointed out before (see Subsection 16. ip the classical approach is not 
convenient. 

We also display an unsteady simulation on FiG. [T71 where the results of the two 
methods are very well reproduced. 



6.3 Remarks concerning mixed still water steady state 

The previous sections deals with the well-balanced property of the numerical 
scheme (|4HI42p -(l44 | for free surface and pressurized flows. To use the well- 
balanced scheme developed in Subsections (|6.ip and (|6.2p . we start from the 
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discrete representation of a mixed still water steady state (as displayed on FiG. 
[TT|) . Assume that there exists io such that all cells rrii on the left hand side of 
the interface are free surface while the other are pressurized: 

( E., = if is^io 
\ E^ = l if i>io ' 

The interface is such that the speed of propagation of the interface is 

nuU: w" = ^ = (see FiG. \m. Therefore, UM = and 

10 + 1/2 An An ^ ' ' ' 

UP = u+. 



U" = UM 


U+ = UP 




c 






FS \ 

/ \ 


/ Press. 

/ \ . 



Figure 13: Mixed free surface-pressurized still water steady state 

For example on FiG. [131 we apply the free surface numerical scheme on the left 
hand side of the interface and the pressurized one on the right cells . As the 
EWBS preserves the pressurized and free surface still water steady state, it also 
preserves the mixed still water steady state. 

7 Numerical tests 

The numerical validation for pipes with constant section and slope has been 
previously studied by two of the authors in [71 [5] and thus are not presented 
in this paper. Since experimental data for mixed flows in any pipe are not 
available, we focus on the behavior of our method for several circular cross- 
section contracting and expanding pipe. Notice that, the equivalent pipe method 
is not relevant for the mixed flows as pointed out by [H [551 130! for instance. 
The mixed flow case is numerically performed on a water hammer test. Start- 
ing from an horizontal free surface still water steady state, the water hammer 
occurs immediately after the increase of the upstream piezometric head while 
the downstream discharge is set to 0. The prescribed hydrograph produces a 
travelling wave which produces a pressurized state propagating from upstream 
to downstream end. Physically an trapped air pocket may appear: it is not 
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taken into account in the PFS-model. Actually, the trapped air pockets van- 
ish or move; some parts of these pockets undergo condensation/vaporisation 
and others parts move and lead to a two phase flow. Consequently the sound 
speed decreases. As our model does not take into account these phenomena, 
the value of c is assumed to be constant. Moreover we should have to deal with 
the entrapment of air bubbles which have a non negligible effect (see [121 for 
instance) . 

The numerical experiments are performed in the case of a 100 m long closed 
circular pipe at altitude bo — Im with slope which corresponds to the el- 
evation and slope of the main pipe axis (we have Z = b{X) = 0, VX). The 
Manning roughness coefficient is l/K^ = 0.012 s/m}-^^. The simulation starts 
from a steady state as a free surface flow with a discharge Q = m? /s. The 
upstream boundary condition is a prescribed hydrograph (see FiG. I19[) while the 
downstream discharge is kept constant to Om'^/s (as displayed on FiG. [19]). We 
compare then the results obtained for uniform, contracting an expanding pipes. 
For each test, the parameters are the same except the downstream diameter: 
the upstream diameter is kept constant to D = Im. The contracting pipe is 
chosen for D — 0.6 m and the expanding one for 13 = 1.4 m (where D denotes 
the downstream diameter). Let us recall that the zero water level corresponds 
to the main pipe axis. The piezometric head is defined as z + p: 



Results are then represented on FiG. [20l The sudden raise of the upstream 
piezometric level produces a pressurized state with a travelling wave. A water 
hammer is then observed since the downstream discharge is null. A careful 

analysis of the flow which is performed by the variable E or equivalently by — 

(see Remark 15.41 and FiG. [TO|l shows that after this transition point, the flow is 
pressurized but in depression which starts approximativcly at time 19s for the 
contracting pipe, 24s for the uniform pipe and 28s for the expanding one. We 
display the piezometric line for different times around the depression time (see 

Fig. [2T|) and also the graph of the function — which confirm that the observed 

times correspond exactly to a depression state for each pipes (see FiG. [Hj). 
We also observe a little smoothing effect and absorption due to the first order 
discretisation type. 

8 Conclusion 

We have derived a free surface and a pressurized model which have been coupled 
using a common set of variables and a suitable pressure law. We have thus 
obtained a mathematical model for unsteady mixed flows in non uniform water 
pipes, that we have called PFS-model. This model takes into account the local 
perturbation of the section and of the slope. 
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The PFS model is numerically solved by a VFRoe scheme using the inter- 
facial upwind to include the source terms into the numerical fluxes. We have 
shown that the classical approximation of the convection matrix (the Roe matrix 
approximation) is not suitable to preserve the still water steady state (except 
for the pressurized case where the value of c helps the scheme to maintain this 
state). Moreover, we have proposed a manner to obtain an exactly well-balanced 
scheme. 

As mentioned in [7] this numerical method with the classical approximation 
of the convection matrix, for constant section, reproduces correctly laboratory 
tests for uniform pipes and can deal with multiple points of transition between 
the two types of flows. As pointed out before, due to the lack of experimental 
data for non uniform pipes, we have only shown the behavior of the piezometric 
line which seems reasonable (at less no major difference was observed). 

As a well-known feature on approximate Godunov scheme, the upwinding of 
the source terms may introduce stationary waves with a vanishing denominator 
when critical flows occurs i.e. u ~ c. Moreover, in its actual form, the presented 
numerical scheme is not able to deal with drying and flooding area. Nevertheless, 
it may be possible to introduce a cut-off function to avoid division by zero 
for each problems: critical stationary waves, drying area and flooding area. 
But it is not the better choice that we can propose, since, truncation of the 
wet area induces a loss of water mass leading to the non-conservativity of the 
mass. Nevertheless, at the present time, we are interested in a mathematical 
kinetic formulation of the PFS model and the construction of a numerical kinetic 
scheme that avoids all these drawbacks |^ . 

The next step is to take into account the air entrainmcnt which may have non 
negligible effects on the behaviour of the piezometric head. A first approach has 
been derived in the case of perfect fluid and perfect gas modelised as a bilayer 
model based on the PFS-model '5^. 

Aknowlegments: The authors wish to thank the two referees for their careful 
reading of the first version of the article and useful remarks. 
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Discharge at x = 40 m 
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Figure 14: The numerical scheme (|4m42|) -(j44 |) with (gS]) and the EWBS for 
pressurized still water steady state with c = 30. 



41 



Discharge at x = 40 m 
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Figure 15: The numerical scheme (|4m42p -(|44 |) with (gB]) and the EWBS for 
pressurized stiU water steady state with c = 200. 
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Discharge at x = 40 m 
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Figure 16: The numerical scheme (|41ll42j) - (|44)) with ([46]) and the EWBS for free 
surface still water steady state. 
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Figure 17: A non stationary test to compare the EWBS and the numerical 
scheme (|4T]|42])-(j44l) with gg. 
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Initial steady state for D = 0.6 m 
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Figure 18: Initial still water steady state for contracting, uniform and expanding 
pipes . 
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Figure 19: Boundary conditions. 
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Figure 20: Piezometric head and discharge a.t X — 50m. 
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Figure 21: Observation of the depression localised approximatively at time t = 
19.117 (contracting pipe), t = 24.075 (uniform pipe) and t — 28.395 (expanding 
pipe). 
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Figure 22: Observation of the depression localised approximatively at time t = 
19.117 (contracting pipe), t = 24.075 (uniform pipe) and t — 28.395 (expanding 
pipe). 



Depression time t = 20. 17 with D = 0.6 m 
9 I 1 1 1 




1 

I 

20 40 60 80 100 

Length (m) 



